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CN \ Abstract: We present results of numerical simulations of the 2-\-ld Nambu - Jona- 
' Lasinio model with non-zero baryon chemical potential n and spatially- varying complex 
Q ■ diquark source strength j. By choosing arg(j) to vary smoothly through 27r across the 

^ . spatial extent of the lattice, a baryon number current is induced which in the high den- 
CS| ! sity phase remains non- vanishing as |j| — >■ 0; we are hence able to extract a quantity 
characteristic of a superfluid known as the helicity modulus. We also study supercur- 
^ ' rent flow at non-zero temperature and estimate the critical temperature at which the 
^ ■ normal phase is restored, which is consistent with the conventional picture for thin-film 
superfiuids in which the transition is viewed in terms of vortex - anti-vortex unbinding. 

^ : 1 Introduction 

^ ' There are unfortunately rather few quantum field theories amenable to study using 



lattice Monte Carlo techniques in the presence of a non-zero chemical potential /i, or more 



O . specifically with ^/T ^ 1. Many important theories, including QCD, cannot be studied 
^ ' because their path integral measure with /i 7^ is not real on analytic continuation to 
Euclidean metric, making Monte Carlo importance sampling inoperative. Of those with 
^ ' positive definite measure, the Nambu Jona-Lasinio (NJL) model with Nf = 2 quark 
fiavors P is one of the most interesting. At /i = the theory exhibits dynamical chiral 
symmetry breaking, with the generation of a constituent quark mass scale S much larger 
than the bare mass m. For > /i^ ~ S, chiral symmetry is restored, and the ground 
state is a degenerate fermi system with fi = Ep c:^ kp [2^. In d+1 dimensions the baryon 
density in this case is ub = ANfi/6{ix — /ic)/(47r)2dr(|). 

The precise nature of the ground state at high density depends on d. For the real- 
istic case c? = 3, lattice simulations suggest that condensation of diquark pairs at the 
Fermi surface takes place leading to spontaneous breakdown of the U(1)b baryon num- 
ber symmetry jH]. An energy gap A > to excite fermionic quasiparticles develops; for 
phenomenologically-motivated lattice parameters the simulations predict A/S ~ 0.15, 
in good agreement with self-consistent model calculations of the gap in superconducting 
quark matter jl]. In this case the NJL model appears to behave as an orthodox BCS su- 
perfiuid; there is long-range ordering of the ground state signalled by the non-vanishing 
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condensate {qq) ^ 0, and a dynamically-generated mass scale A. Since a U(1)b sym- 
metry has been spontaneously broken, there is a massless scalar qq bound state in the 
spectrum, which is associated with long-range interactions between vortex excitations 
in the superfluid, and with a collective propagating mode for T > known as second 
sound. 

However, both for obvious numerical convenience, and for a more formal reason, 
namely the existence of an interacting continuum limit, lattice simulations were first 
applied to the NJL model with 7^ in 2-1-1 dimensions IHl El- In this case the 
physics appears radically different. Whilst there is evidence for long-range coherence 
of diquark correlation functions [Sj, there is no long-range order, and apparently no 
gap. Rather, the condensate vanishes non-analytically with the diquark source strength, 
{qq) oc with < a{jji) < 1 0- The results were interpreted in terms of a critical phase 
for all n > fic, in which the diquark correlator decays algebraically, {qq{0)qq{r)) oc 
[7j. Since all simulations are performed on finite systems, and therefore necessarily at 
T > 0, the absence of long-range order is consistent with the Coleman- Mermin- Wagner 
theorem for 2d systems 0. The situation is analogous to the low-T phase of the 2d 
X-Y model, one of whose physical applications is the description of superfluidity in thin 
films P?. 

The defining property of a superfluid is that the flux density of conserved charge, or 
supercurrent J, is related to the spatial variation of the phase angle 9 of the U(l)-valued 
order parameter field (in this case (qq)) via 

J = TW. (1) 

The constant of proportionality T is known as the helicity modulus. For a textbook 
non-relativistic superfluid such as ^He it is given by 

= (2) 

where M is the mass of the helium atom and Ug is a parameter called the superfluid 
density, which need not coincide with the charge density of the condensate. For a 
relativistic system T is best thought of as a phenomenological parameter in its own 
right, rather like in {d + l)-dimensional chiral model ^U]. One way of understanding 
superfluidity in the absence of long-range order in a 2d system is to observe that the 
only way to change the quantised circulation § J.dl around one direction of a finite torus 
is to excite a vortex - anti-vortex pair, and transport one member of the pair around 
the other direction of the torus before re-annihiliation. The energy required to do this 
scales as InL where L is the size of the system j^]: hence in the thermodynamic limit 
circulation patterns are topologically stable. 

This Letter will present further support for this scenario in NJL2+1 by extracting 
T via a calculation of the induced baryon number current J = (i/jjip) in response to a 
spatially varying diquark source. As well as providing direct verification of superfluid 
behaviour in a fermionic model, we will also study the behaviour of J as temperature 
T is increased, and find the transition to "normal" behaviour at a critical Tc consistent 
with analytic expectations. 
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2 Method 



The lattice NJL model studied is identical to that of [3111: 

Snjl = E XxM^^yXy + JX. ^2Xx + JX^raX. + 4 E t^^l** > (3) 

X 9 X 

where Xi X ^"^^ isodoublet staggered lattice fermion fields, $ = (xl +i7f.r is an auxiliary 
bosonic field defined on the dual lattice sites x, and the matrix M is 

v=0,l,2 ^ I <x,x> 

(4) 

Here < x,x > denotes the set of 8 dual sites x surrounding x, ri^^ = jg 
the Kawamoto-Smit phase required for a Lorentz covariant continuum limit, and Ex = 
^_lY'0+xi+x2 _ A full description of the symmetries of (jS)) and the numerical simulation 
method is given in [7]. The only novelty in the current study is that the diquark source 
strengths j, J are now specified to be spatially varying, or "twisted": 

j = jo exp(i6'^) ; j = jo exp(-i6'^) (5) 

with jo a real constant. To ensure homogeneity and single- valuedness on an x Lt 
lattice we demand 

271 ^ 27r 

9 = —{riiXi + n2X2) =^ V6' = — (^i, na). (6) 

A constant supercurrent of the form (Q) is therefore specified by a pair of integers (ni, na). 
It remains to define the conserved baryon number current J^: 

Jux = lie^'-'^XxXx+p + e-'^'-^'XxXx-o). (7) 

The timelike component of (|7j) is none other than the baryon charge density hb reported 
in [3[7]. Here we shall use the same stochastic technique to estimate the quantum expec- 
tation value of the spacelike components /i) = (L^Lt)'^ Jx{j, fi) to demonstrate 
behaviour of the form The strategy will be to compute J for fixed (^1,77,2) for a 
range of jo and extrapolate Jo ^ 0. Behaviour consistent with in this limit is deemed 
to be superfluid. 

We used the same simulation parameters as jSJ [7], namely = 2.0a, ma = 0.01, 
which at /i = yields a dynamically-generated constituent mass, which in effect sets the 
scale, of Sa = 0.71. As /i is raised, there is a sharp transition from a chirally broken 
vacuum with {xx) — '^^1 ub — to a. chirally restored phase with > at ~ 0.65. 
Studies of the fermion dispersion relation in the phase ii> Hc are consistent with a sharp 
Fermi surface with kp^fJ' and vanishing gap A ~ (HllZl- 
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3 Results 
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Figure 1: J2 vs. jo on a 16^ x Lt lattice for two different values of fi. 

In Fig. ^ we plot J2 (strictly its imaginary part) as a function of jo for lattices of 
various temporal extent Lf at two representative values of /i: fia = 0.2 lies in the chirally- 
broken low density phase, and /la = 0.8 in the high density phase, where n^a^ — 0.25 
In all the plots shown here we have chosen (ni,n2) = (0,1) to minimise lattice 
artifacts, and from now on we set the lattice spacing a = 1. 

The contrast between the two phases is quite dramatic. For yU = 0.2 J2 appears to 
vary approximately quadratically with jo, and extrapolate to zero as jo 0. There is 
no significant effect as Lt —>■ 00, or alternatively as T ^ 0. At = 0.8 the small- jo 
behaviour depends very sensitively on Lf, as T ^ the data accumulate on a straight 
line which clearly extrapolates to a non-zero value as jo — >■ 0. 

This behaviour is readily explained by writing the order parameter (diquark) field 
as = 00^*^, with 00 approximately constant. A natural effective Hamiltonian for long 
wavelength order parameter fluctuations at low temperature is then 



eff 



-(V0)*.(V0) ~ ^{V9f 



(8) 



The corresponding Noether current is J 



^V0 - (V0* 



blVO. For ^ < fic, 
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it is natural to postulate proportional to j, leading to J2(jo 



assume that lim 



7^ we recover with T 



oc Jq. For /X > /ic, if we 
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Figure 2: J2 vs. jo on a 16^ x 64 lattice for various fi. 

With confidence that = 64 suffices to determine the phase, in Fig. |21we plot J2U0) 
for various fi. There is a sharp change between values of < 0.65, which display the 
low density quasi-quadratic behaviour and smoothly extrapolate to zero as jo 0, and 
fi > 0.68 which show a negative curvature characteristic of the high density phase. We 
thus determine the critical chemical potential for the onset of superfluidity 0.65 < fic < 
0.68, in good agreement with the critical value for chiral symmetry restoration. Since 
as yet we have no systematic method of extrapolating to jo — > for /x^/Xc to obtain an 
estimate for J2{fJ') as an "order parameter", we can make no decisive statement about 
the nature of the transition, but note that the behaviour of J2(jo) varies much more 
sharply across the transition than the diquark condensate (gg+(j)), either in this model 
(Cf. Figure 2 of [3), or even in NJL3+1 (Cf. Figure 4 of [3J). This matches the sharp 
drop in the chiral order parameter {xx) and corresponding rise in at /i = /ic j3IZ|, 
and is consistent with the analytic prediction of a strong first order transition in the 
large- A^/ limit [TT]. 

From now on we work exclusively at fi = 0.8 in an attempt to understand the 
superfiuid phase further. In Fig.Elwe plot T = J2-^s/2vr versus jo on x64 lattices. Just 
as for the {qq+{j)) data |7|, it turns out that the data are well-fitted by T(jo) = A+B/Lg, 
resulting in the Ls ^ 00 extrapolation shown in the plot. Recall that as well as genuine 
finite-size effects in this case there may also be some discretisation artifacts, since as 
Lg increases the gradient operator in (jH)) becomes better-approximated by the finite 
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Figure 3: T vs. jo on a x 64 lattice for /i = 0.8. 



difference. Finally the T(jo) data in the thermodynamic limit are extrapolated to jo 
with a remarkably simple linear fit, resulting in T = 0.1413(14). We thus quote a result 
for the hehcity modulus of T/S = 0.200(2) at /xa = 0.8. 

It is interesting to pause and ask what value might be expected for T in a conventional 
symmetry- breaking scenario. Let us define diquark operators qq± = \{x^^^2X i X'^2X^^) 
and source strengths j± = j ^ J, so that the diquark terms in the action Q read 
~^ j-QQ-- III thfi limit j_ = the equation of motion for the current is then 

A-J^ = 2uqq_. (9) 

In the same limit the U(l)B-equivalent of the axial Ward identity reads 

= j+E(9?40)gg-(^)) = j^mqq-Hl'- (lo) 

where the second equality assumes that the correlation function is dominated by a 
pseudo-Goldstone pole of the form (/c^ + M^)"^, and |— ) denotes a one-Goldstone state. 
We now introduce the U(l)B-equivalent form of the PCAC hypothesis: 

m-J,\-) = v/TAf! = 2jV(0|gg-|-) (11) 

where we have used the relation T = derived in ^U] , and the second equality follows 
from Q. Combining (|10|) and (|11|) leads to the equivalent of the "Gell-Mann-Oakes- 
Renner" relation: 

TGMORM^_ = 8j{qq+). (12) 
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This can be compared with numerical data for {qq+{j)) and M^{j) in [7j. At j = 0.3, 
M_ = 0.95, {qq+) = 0.72 yielding Tgmor ~ 1-9; at j = 0.1, M„ = 0.4, = 0.52 

yielding Tgmor ~ 2.6. We conclude T <^ Tgmor, consistent with our hypothesis 
that no symmetry breaking occurs, but that the dynamics are dominated by long-range 
phase fluctuations of the order parameter field, described by a strongly-interacting scalar 
diquark field rather than a weakly-interacting Goldstone mode. 

3.2 T > 

In this section for the first time we explore the superfluid phase at non-zero temperature. 
We expect a restoration to the normal phase at some critical T^. In a comparable 
numerical study of NJLa+i which exhibits superfluidity via orthodox symmetry breaking 
[2], the value of could be estimated from the zero temperature gap A via the BCS 
prediction A/Tc ~ 1.76. Since this implied that Lt had to exceed 35a for the system to be 
superfluid, an unambiguous extrapolation j ^ to permit a systematic study of T > 
was not possible. In the current case we shall see that although the j — > extrapolation 
still remains a problem, attaining T < is well within reach of the simulation. 

First let us review a heuristic argument for the expected value of T^, starting from 
the Hamiltonian Hgff (jHl) with (pl = T The phase field 9{x) may be disrupted by 
topologically non-trivial vortex excitations of the form 6 = qip, \V6\ = q/r, where q is 
integer and x is written {r,il)). The energy of a single vortex is thus 

^~I^^'27rrc/r(^^)' = 7rTg2ln(^:^^. (13) 
Since a vortex can be located on any of lattice sites, the entropy 

S = 2\n(^^y (14) 
The free energy F = E — TS thus changes sign for q = I vortices at a critical temperature 

Tc = ^T. (15) 

The interpretation is that a phase transition separates a low-T superfluid phase in which 
vortices are confined to bound dipole pairs, and a high-T normal phase in which the vor- 
tex anti- vortex plasma screens the long-range interactions responsible for the divergent 
energy in (jl3p . The relation (jl5|) remains valid in a more sophisticated renormalisation 
group treatment, except that T is now T-dependent and should be replaced by its value 
T(Tc) exactly at the transition [T^ . 

Combining our result for T with ()15|) yields an estimate Lf ~ 4.5 for the temporal 
lattice extent where a transition to the normal phase might be expected at /i = 0.8. Fig.|3] 
shows J2(jo) on 32^ x Lf lattices with Lt ranging from 64 all the way down to 2. At the 
extremes Lf > 32, Lf < 4 the data are reminiscent of those characterising respectively the 
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Figure 4: J2 vs. jo on a 32^ x Lf lattice at = 0.8 for various Lt. 

high and low baryon density phases in Fig. ^ For intermediate temperatures, however, 
J2U0) shows positive curvature near the origin followed by negative curvature at larger 
jo, and once again the means of extrapolating jo — > to determine whether superfluidity 
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Figure 5: J2 vs. Lt on a 32^ x lattice at /x = 0.8 for fixed jo = 0.025. 



persists is unclear. 

In Fig. El we try a different tactic, plotting J2 for every even G [2,64] for fixed 
jo = 0.025. A linear fit Jo = a^Lt + ai through data with Lt < 42 seems quite reasonable, 
yielding ao = 0.00212(25), oi = —0.01537(9). If we identify the intercept on the Lt- 
axis with the transition, we deduce Ltc = 7.25(95) and hence T/Tc = 1.02(13), to be 
compared with the theoretical value 0.637 from |T3|). 

4 Summary 

In this short study of the response of the system to a twisted diquark source forcing a 
baryon number current, we have provided direct evidence for the superfiuid nature of the 
ground state of NJL2+1 at high baryon density, and quantified it at one representative 
value of /i via the helicity modulus T. It should be stressed that the "physical" value 
T /Ti 0.2 quoted is still to be extrapolated to the continuum limit. It is probably 
more important to note that the numerical value of T is an order of magnitude smaller 
than might be expected in an orthodox symmetry breaking scenario, and is consistent 
with the non-Goldstone, strongly self-interacting nature of the scalar diquark excitations 
above the ground state. 

We also studied the response of the system to non-zero temperature. Whilst we were 
unable to extrapolate to the zero source limit in a controlled way, by studying fixed 
jo we were able to estimate a critical temperature Tc for breakdown of superfluidity of 
the same order of magnitude as, and only slightly smaller than, the Kosterlitz-Thouless 
prediction for a 2d system with U(l) global symmetry, which follows from characterising 
the superfluid/normal transition as arising from vortex pair unbinding. More refined 
simulations would be required to determine whether actually has the KT value, or 
whether NJL2+1, which in addition to the scalar diquark excitations contains massless 
fermion degrees of freedom, actually hes in a different universality class, as suggested by 
estimates of the critical exponent S |7] . 
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